Bose-Einstein condensate collapse: a comparison between theory and experiment 
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We solve the Gross-Pitaevskii equation numerically for the collapse induced by a switch from pos- 
itive to negative scattering lengths. We compare our results with experiments performed at JILA 
with Bose-Einstein condensates of 85 Rb, in which the scattering length was controlled using a Fesh- 
bach resonance. Building on previous theoretical work we identify quantitative differences between 
the predictions of mean-field theory and the results of the experiments. Besides the previously re- 
ported difference between the predicted and observed critical atom number for collapse, we also find 
that the predicted collapse times systematically exceed those observed experimentally. Quantum 
field effects, such as fragmentation, that might account for these discrepancies are discussed. 

PACS numbers: PACS numbers: 03.75.Fi, 03.75.Be 



Introduction - Most experiments on dilute gas Bose- 
Einstein condensates (BECs) are performed with atoms 
that have a repulsive two-body interaction. Exceptions 
are the experiments on 7 Li |l], || and, more recently, on 
85 Rb ||, ||]. For 85 Rb a Feshbach resonance allows the 
two-body interaction strength to be tuned over a wide 
range of attractive and repulsive values. In particular, 
the scattering length has been rapidly switched from pos- 
itive (repulsive interaction) to negative (attractive inter- 
action) values, leading to the collapse and subsequent 
explosion of the condensate. Recently, the large positive 
scattering lengths attainable in this system have been 
used to produce atom-molecule condensates ||. 

In the following we report on our modelling of the 85 Rb 
collapse experiments ||, using the Gross-Pitaevskii (GP) 
equation for the expectation value of the field operator 
§ @, §, §. Saito and Ueda fo| and Adhikari gj have 
also modelled these experiments by numerical solution 
of the cylindrically symmetric GP equation. Saito and 
Ueda conclude that this describes the collapsing and ex- 
ploding dynamics at least qualitatively jIo|. Following 
their suggestion, we report a more quantitative compari- 
son between the theoretical and experimental results, and 
find significant differences. 

The series of experiments on the collapse and explosion 
of 85 Rb BECs challenges theoretical models in a number 
of ways iQ. A body of theoretical work based on the 
GP equation predicts the critical number of atoms N cr 
for collapse to be significantly larger than is observed. 
The expression for the critical number is 
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where a^o = \fh/(rnuj) is the harmonic oscillator scale 
length, with u> the geometric mean of the trap frequencies 
in the three Cartesian directions, and a the scattering 
length. Experimentally, k — 0.46 ± 0.06 [0], whereas 
k = 0.57 according to various approximate solutions of 
the GP equation (r| 

We have confirmed this GP prediction for the specific 



cylindrically symmetric experimental case [|J with cylin- 
drically symmetric numerical solutions. We verified these 
with full three dimensional numerical solutions, and also 
confirmed that slight departures from cylindrical sym- 
metry had no effect on the critical number p5[ . Con- 
sequently there is a disagreement at the two standard 
deviations level, which should be regarded as significant. 

We also report a new quantitative discrepancy between 
the predictions of the GP model and experiment. Under 
certain conditions, the GP predicted time to the initi- 
ation of collapse, ^collapse, is systematically longer than 
that observed in the experiments Q. 

The GP model - In the conclusion we will discuss the 
possiblility that these discrepancies result from quantum 
field effects beyond the GP approximation. We there- 
fore now derive the GP equation from the quantum field 
theory. 

The second- quantised Hamiltonian for a dilute gas, in 
terms of the field operator fy(r,t), is 



H = J dr&H Q y 

+ ~ J drdr'&&'V(r - r')*'$, (2) 

where ifr' = ^>(r',t) and Hq is the single particle Hamil- 
tonian for the kinetic energy and trapping potential 
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where m is the atomic mass (1.41 x 10 -25 kg for 85 Rb), 
and u>i is the trap frequency along Cartesian axis i. In 
the limit of particles separated by distances much greater 
than the scattering length a we approximate the two- 
body potential by a delta function interaction [§[ 0, |§[ ||] 



V(r-r')=gS(r-r'), g 
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The Heisenberg dynamical equation for the field operator 
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is then 



ih—J> = Hq^ 
at 



(5) 



Taking the symmetry-breaking approach we assume that 
the field expectation value is not zero and define it as the 
GP wavefunction (\&(r,i)) = <I>(r,i), normalised to the 
number of particles N 



N 



|$(r,t)| 2 dr. 
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Then taking the expectation value of the Heisenberg 
equation (||) gives 



d 
ot 
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If we assume that the expectation value factorises, as it 
would, for example, if the system were in an eigenstate 
of the field operator, 



then we obtain the GP equation 

ih^ = (H +g\<S>\ 2 )<S>. 
ot 
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In order to model atom loss due to three-body recombi- 
nation we add a phenomenological term proportional to 
the density |$| 2 squared with rate coefficient if 3/2 M llq| 



ih^-$=(H +g\<!>\ 2 
at 
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We assume one-body and two-body loss are negligible, 
as was true for the relevant experiments. The number of 
atoms then decays as 



^ = -K 3 J |$M)| 6 dr. 



(11) 



GP Results - As an example of the ability of the GP 
equation to correctly model the 85 Rb || experiments we 
present Fig. [l]. It is the result of a numerical solution of 
the (two dimensional) cylindrically symmetric GP equa- 
tion for <l(r, z) 



h 2 

+ -m(uj 2 r 2 + w 2 z 2 )$ -I- g|$| 2 $ 
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Parameters are the same those of Fig. lb of Donley et 
al. ||. Specifically, the ground state of the GP equa- 
tion for a = +7ao was switched in 1 ms to a = — 30ao, 
where ao = 0.0529 nm is the Bohr radius. For the three- 
body recombination rate coefficient K3 — 190 x 10~ 28 
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FIG. 1: Experimental and numerical results for the number 
of atoms N versus time after a switch from a — +7ao to 
a — — 30ao- The experimental points (x) are from Fig. lb of 
Donley et al. Ji). The numerical results are for Kz — 190 x 
10~ 28 cm 



(filled circles) and for K 3 = 78 x 10" 



+ 



Other parameters are as given in the experimental paper 
iVo = 16, 000, radial frequency ui r = 2tv x 17.5 Hz, axial 
frequency uj z — 2ir x 6.8 Hz. 



cm 6 s _1 the agreement with the experimental results is 
good. However it should be noted that the experimen- 
tal points are the "remnant" atom number, while the 
numerical points are the total atom number, which over- 
estimates the remnant atom number. A smaller value of 
-K3 agrees better with the earlier points, while overesti- 
mating the final atom number. The precise value of K3 
has little effect on the conclusions of this paper, which 
concern the initiation of collapse. 

These results agree with those reported by Saito and 
Ueda 10 and Adhikari jll[] . However the former authors 
used a much smaller value of the three-body recombina- 
tion rate coefficient K3 = 2 x 10~ 28 cm 6 s _1 . This pro- 
duces the collapses and revivals in condensate size that 
were observed in their simulations. These only become 
important for K3 less than about 10 -26 cm 6 s _1 . Ad- 
hikari |llj used the much larger value if 3 = 13 x 10~ 25 
cm 6 s -1 . Since three-body recombination is responsible 
for the atom loss, it is remarkable that such a wide range 
of coefficients reproduces the experimental results. 

The three-body recombination rate coefficient K3 is 
expected to vary strongly near the Feshbach resonance 
J17J. Experimental determination of K3 is difficult due 
to the low densities of 85 Rb condensates. Upper bounds 
have been estimated to be 5 x 10 -25 cm 6 s _1 , dropping 
to 10~ 26 cm 6 s _1 nearer the Feshbach resonance Jlq]. 

The cylindrically symmetric numerical simulations 
were performed on a 512 x 512 grid, 35.64 /mi long in 
the axial (z) direction and with the radial coordinate 
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extending to 11.88 /im. The corresponding spatial grid 
spacings were therefore 0.07 /jm and 0.023 fim. The time 
steps were 2.34 ns. All simulations were performed on 
a multiprocessor machine ]l8[ ], using up to 32 proces- 
sors, and the RK4IP algorithm developed by the BEC 
theory group of R. Ballagh at the University of Otago 
p9[ . This is a pseudo-spectral method with a Runge- 
Kutta time step. The cylindrically symmetric and full 
three-dimensional codes are independent and were cross 
checked. Grid spacings and time steps were varied to en- 
sure convergence. Overall the results were found to be 
quite robust. As another test, we solved the GP equation 
for a half a radial period after the quenching of the col- 
lapse. As was observed experimentally, the condensate 
refocussed onto the axis, due to the oscillation in the 
harmonic radial potential. All this, together with the 
agreement of our results with those of Saito and Ueda 
fiof and Adhikari 11 , gives us confidence in their ac- 
curacy. Following Adhikari pj| , the initial condition for 
Fig. |l| was generated by adiabatically expanding the har- 
monic oscillator initial state a = to a = +7ao over 444 
ms. 

Figure || presents our calculations of the collapse times 
^collapse for the conditions of Fig. 2 of Donley et al. ||. 
The collapse times were determined by visually fitting 
plots of atom number versus time to the functional form 

N = (No - Nf) exp[-(i - t C ollapsc)/T d ccay] + Nf, (13) 

where Nf is the long time atom number. An example 
is given in the inset to Fig. |^. We have also plotted 
the experimental results reported in Fig. 2 of Donley et 
al. ||, and find a small, but significant, systematic dis- 
agreement with the GP results. Although the reported 
errors in the experimental collapse times are large, the 
GP values for i C oiiapse are consistently longer than the 
experimental ones. This is surprising as the GP model is 
expected to be valid for the low densities preceding the 
collapse. If it were to fail, it would be expected to do so 
at the high densities generated subsequently. Neverthe- 
less, the disagreement is not unprecedented since, as we 
discussed earlier, the GP model also overestimates the 
critical number for collapse. 

The estimates of i C oiiapse by Saito and Ueda |l(J (their 
Fig. 3) are between five percent (low a) and ten per- 
cent (high a) smaller than ours. This is consistent with 
the smaller three-body recombination rate coefficient K3 
they used. However their results are still significantly 
longer than the experimentally measured times. 

We have confirmed these cylindrically symmetric simu- 
lations by performing full three-dimensional simulations. 
In particular we broke the cylindrical symmetry by using 
trap frequencies of 17.24 x 17.47 x 6.80 Hz j|. 

We were unable to substantially improve the agree- 
ment by either changing the initial condition to reflect 
the experimental uncertainty of a = ±2ao, or by varying 
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FIG. 2: Experimental and numerical results for the collapse 
time Collapse versus scattering length a co n apsG after a switch 
from a = to a co ii apsc . The experimental points £+) and 
their error bars are from Fig. 2 of Donley et al. ]3|. The 
numerical results (filled circles) are for K3 — 190 x 10~ 28 
cm 6 s _1 . Other parameters are as given in the Fig. [l] caption, 
except: iVo = 6,000. Inset: Example of the fitting procedure 
used to determine the collapse times. Shown is a fit of the 
functional form Eq. ( p"3| ) (solid line) to the GP simulation ( x ) 
for a — — lOao. The fit parameters here are t co ii apso = 9.8 ms, 
Tdocay = 0.7 ms, and Ni/No = 0.5544. 



the three-body recombination rate coefficient. This sug- 
gests that some of the physics determining the collapse 
time is not captured by our GP model. 

Discussion.- Both the collapse time and critical num- 
ber discrepancies could be resolved by using a scatter- 
ing length in the GP model larger in magnitude than 
the experimental value. This would reduce the collapse 
time and decrease the critical number, as required. The 
required increases in the scattering length magnitudes 
vary, ranging from a factor of 0.57/0.46 — 1.2 for the 
critical number, up to a factor of about two for the col- 
lapse times for large a co iiapse- However, the scattering 
length is experimentally well calibrated p0[ , so any such 
change would reflect a deficiency of our GP model. 

One possible origin of the discrepancy is the effect 
of thermal non-condensed atoms. Because of the quan- 
tum statistics of collisions betwen bosons, the scattering 
length between a condensed atom and an atom in another 
mode is twice that between two condensed atoms. Hence 
one might expect the presence of thermal uncondensed 
atoms to shorten the collapse time compared to the GP 
prediction, as observed. Furthermore, this might be ap- 
proximately corrected for by using an increased magni- 
tude effective scattering length in the GP model. How- 
ever the uncondensed fraction is much less than 10% of 
the total number of atoms so it seems unlikely that 
this effect is large enough to account for the discrepancy. 
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Furthermore, Roberts et al. |5| reported that the critical 
number for collapse N C1 - was insensitive to varying the 
temperature. Therefore we do not expect finite tempera- 
ture extensions of the GP theory to explain the discrep- 
ancy f§. 

Another possibile origin is quantised atom field effects. 
These might arise due a breakdown of the factorisation 
assumption Eq.(||). There have been several suggestions 
for how the quantised field might influence the collapse 
p2| , |23| , p4j . Furthermore, Nozieres J25| has emphasised 
that only for positive scattering lengths does an energy 
barrier protect BECs from fragmentation into many pop- 
ulated states. For negative scattering lengths, mean-field 
energy is released when atoms scatter from the conden- 
sate into other modes. Fragmentation could increase the 
effective scattering length by up to a factor of two. 

In order to investigate the behaviour of a fully quan- 
tised atom field, we have used the gauge-P function ap- 
proach recently developed by Deuar and Drummond [p6| . 
This method overcomes some of the problems that plague 
stochastic simulations based on the positive P-function 
quasi-probability distribution p7[ |28fl . We were compu- 
tationally limited to simulations in one spatial dimension 
and found agreement with the GP collapse times at the 
one percent level. 

Although this preliminary work does not provide evi- 
dence for quantum field effects, it is important to extend 
the fully quantised field modelling to three spatial dimen- 
sions, and hence to use actual experimental parameters. 
As shown in Fig. 1, there are parameters for which the 
GP theory does agree with experiments. One approach 
is the recently developed perturbation theory which ex- 
tends the GP model to include normal and anomalous 
densities of the quantum field |^] . This method has re- 
cently been successfully applied to the formation of atom- 
molecule condensates in 85 Rb fl30| . 

This research was supported by the Australian Re- 
search Council and by an award under the Merit Allo- 
cation Scheme on the National Facility of the Australian 
Partnership for Advanced Computing |18] . 
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